Direct correlation functions of the Widom-Rowlinson model 
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We calculate, through Monte Carlo numerical simulations, the partial total and direct 
correlation functions of the three dimensional symmetric Widom-Rowlinson mixture. 
We find that the differences between the partial direct correlation functions from 
simulation and from the Percus-Yevick approximation (calculated analytically by 
Ahn and Lebowitz) are well fitted by Gaussians. We provide an analytical expression 
for the fit parameters as function of the density. We also present Monte Carlo 
simulation data for the direct correlation functions of a couple of non additive hard 
sphere systems to discuss the modification induced by finite like diameters. 
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I. INTRODUCTION 

Fluid binary mixtures may exhibit the phenomenon of phase separation. The simplest 
system able to undergo a demixing phase transition is the model introduced by Widom and 



Rowlinson some years ago ^ . Consider a binary mixture of non-additive hard spheres. This 
is a fluid made of hard spheres of specie 1 of diameter .Rn and number density pi and hard 
spheres of specie 2 of diameter R22 and number density p2, with a pair interaction potential 
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between species i and j that can be written as follows 

oo r < R 



v^Ar)={ ' , (1) 

r > Rij 

where R12 = {Ru + -R22)/2 + a. The Widom-Rowlinson (WR) model is obtained choosing 
the diameters of the spheres equal to 0, 

Rn = R22 = , (2) 

so that there is no interaction between like spheres and there is a hard core repulsion of 
diameter a between unlike spheres. The symmetry of the system induces the symmetry of 



the unlike correlations [hi2{r) = h2i{r), 
in the past by exact i2| and approximate 



r)^^2i{r)]. The WR model has been studied 
methods and it has been shown that it 



exhibits a phase transition at high density. More recently, additional studies have appeared 
and theoretical predictions have been confirmed by Monte Carlo (MC) computer simulations 



i],y,i9 



id 



In this paper we will study the three dimensional symmetric Widom-Rowlinson mixture 
for which pi = p2 = p/2, where p is the total number density of the fluid, and 

hn{r) = h22{r) , (3) 
Cii(r) = C22(r) . (4) 

Moreover we know from (Q) that the partial pair correlation function gij = hij + 1 must obey 

gij{r) =0 for r < Rij . (5) 

Our main goal is to focus on the direct correlation functions (dcf) of the WR model as a 
simplified prototype of non-additive hard spheres (NAHS) systems. The reasons to focus on 
the dcf 's is twofold: on the one hand, they are easier functions to model and fit. On the other 
hand, they play a central role in approximate theories like the Percus-Yevick approximation 
or mean spherical approximation (MSA) We hope that a better understanding of the 
dcf's properties in the WR model, could help in developing accurate analytical theories for 
the NAHS systems. 

We calculate through Monte Carlo simulations the like g^'"'^ (r) and unlike g^'"'^ (r) pair 
distribution functions for a system large enough to allow a meaningful determination of the 
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correspondent partial direct correlation functions c^\^^\r) and c^^'"\r), using the Ornstein 



Zernike equation 



We compare the results for the unlike direct correlation function with 



the results of the Percus-Yevick (PY) analytic solution found by Ahn and Lebowitz 0,0]. In 
the same spirit as the work of Grundke and Henderson for a mixture of additive hard spheres 
[l^ . we propose a fit for the functions ^\i{r) = c[*^'"^(r) and A^gl'^) = ~ c[^^''(r). 

At the end of the paper we also show the results from two Monte Carlo simulations on a 
mixture of non-additive hard spheres with equal diameter spheres Ru = R22 = -R12/2 and 
on one with different diameter spheres Ru = and R22 = Ru to study the effect of non 
zero like diameters on the WR dcf 's. 



II. MONTE CARLO SIMULATION AND PY SOLUTION 

The Monte Carlo simulation was performed with a standard NVT Metropolis algorithm 
[3| using = 4000 particles. Linked lists [l^ have been used to reduce the computational 
cost. We generally used 5.2 x 10^ Monte Carlo steps where one step corresponds to the 
attempt to move a single particle. The typical CPU time for each density was around 20 
hours (runs at higher densities took longer than runs at smaller densities) on a Compaq 
AlphaServer 4100 5/533. 

We run the simulation of WR model at 6 different densities p = pa^ = 0.28748, 0.4, 0.45, 
0.5, 0.575, and 0.65. Notice that the most recent computer simulation calculations P. 
give consistent estimates of the critical density around 0.75. Our data at the highest density 
(0.65) are consistent with a one phase system. 

The Monte Carlo simulation returned the gij{r) over a range not less than 9.175a for 
the densest system. In all the studied cases the pair distribution functions attained their 
asymptotic value well inside the maximum distance they were evaluated. Thus, it has been 
possible to obtain accurate fourier transforms of the correlation functions [hij{k)]. To obtain 
the Cjj(r) we used Ornstein-Zernike equation as follows 

[l + f/^ii(A;)]'- [f/ii2(fc)] 
ci2ik) = ^# , (7) 

[l + Eh,,(k)Y-[lh2ik)Y 

From the hij[k) and Cij{k) we get the difference jijik) = hij{k) — Cij{k) which is the fourier 



4 



transform of a continuous function in real space. So it is safe to transform back in real space 
[to get 7jj(r)]. Finally, the dcf's are obtained from the differences hij{r) — jijir). 

While for a system of non-additive hard sphere in three dimensions a closed form solution 
to the PY approximation is still lacking, Ahn and Lebowitz have found an analytic solution 
of this approximation for the WR model (in one and three dimensions). 

The PY approximation consists of the assumption that Cjj(r) does not extend beyond the 
range of the potential 



Cij[r) 



for r > Rij . (8) 

Combining this with the exact relation ^ and using the Ornstein-Zernike equation we are 
left with a set of equations for Cjj(r) and gij{r) which have been solved analytically by Ahn 
and Lebowitz. 

Their solution is parameterized by a parameter zq. They introduce the following two 
functions of Zq (which can be written in terms of elliptic integrals of the first and third kind) 

dz 



oo 



z^y z^ + Az/zo 
dz 



^Jz^ + ^zjzQ - 4 

and define Zq in terms of the partial densities p\ and as follows 



7] = 2ny/pip2 



cos/i 



(9) 
(10) 

(11) 



They then define the following functions (note that in the last equality of equation (3.76) in 
{4] there is a misprint) 



Cl2(fc) 



X sm 




l + Y 



z^Y^ + 4F + 4 



dz 



{z + Y)^Jzlz^ + Az-A 



huik) = Ci2{k)[l - pip2Ci2{k)] , 



(12) 
(13) 



where Y = {2k/ hf. 

We also realized that some other misprint should be present in the Ahn and Lebowitz 
paper since we have found empirically that the PY solution (with k in units of a) should be 
given by 



Ci2{k) = cuiks) 



(14) 
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where s is a scale parameter to be determined as follows 

s = -[hi2{r = 0)]'/' (15) 

Notice that for the symmetric case pi = P2 = p/2 and rj = np = 0.90316 ... we find Zq = 1 
and s = 1. 

In Figs. andlHlwe show three cases corresponding to the extreme and one intermedi- 
ate density. In the figures, we compare the MC simulation data with the PY solution for the 
partial pair distribution functions and the partial direct correlation functions. Our results 
for the partial pair distribution functions at pa^ = 0.65 are in good agreement with the ones 



of Shew and Yethiraj ^ . The figures show how the like correlation functions obtained from 
the PY approximation are the ones that differ most from the MC simulation data. The 
difference is more marked in a neighborhood of r = and becomes more pronounced as the 
density increases. 



III. FIT OF THE DATA 

From the simulations we found that c^u'"\r) < 8 x 10^^ for r > a at all the densities 
studied. This allows us to say that A^gl'") — r > a. Moreover we found that both 
Al2{r) for r < a, and Al^^r) are very well fitted by Gaussians 

A5i(r) ~ bu exp[-aii(r + rfn)^] for all r > 0, (16) 
Al2ir) ~ 6i2 exp[-ai2r^] for < r < a, (17) 

In Figs. EJandElwe show the behaviors of the parameters of the best fit (fTT)|l and ()17j) . with 
density. In order to check the quality of the fit, we did not use the data at p = 0.45 in the 
determination of the parameters. The points for ai2 and 612 are well fitted by a straight line 
or a parabola. As shown in Fig. |3]the best parabolae are 



ai2(p) = 0.839 + 0.096p - 1.287p2 , 
6i2(p) = -0.155 + 0.759p-0.159p2 



(18) 
(19) 
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Fig. El shows how the parameters for A^]^(r) are much more scattered and hard to fit. The 
quartic polynomial going through the five points, for each coefficient, are 

an(p) = -55.25 + 504.8p- 1659.p2 + 2364.p3- 1236./ , (20) 
6ii(p) = 171.4 - 1556.p + 5166.p2- 7421.p=^ + 3906.p^ , (21) 
duip) = 128.9 - 1144.p + 3747.p2- 5328.p2 + 2782.p^ , (22) 

The difficulty in finding a good fit for these parameters may be twofold: first we are fitting 
A^]^(r) with a three (instead of two) parameters curve and second the partial pair distribution 
functions obtained from the Monte Carlo simulation are less accurate in a neighborhood 
of the origin (due to the reduced statistics there). This inaccuracy is amplified in the 
process of finding the partial direct correlation functions. Such inaccuracy will not affect 
significantly A^2('") which has a derivative very close or equal to zero near the origin, but it 
will significantly affect A^]^(r) which is very steep near the origin. 

In order to estimate the quality of the fit we have used the simulation data at p = 0.45. 
From Fig. |3] we can see how the parabolic fit is a very good one. In Fig. |31 the point at 
p = 0.45 gives an indication of the accuracy of the quartic fit. We have also compared the 
pair distribution and direct correlation functions obtained from the fit with those from MC: 
both the like and unlike distribution functions are well reproduced while there is a visible 
discrepancy in the dcf from the origin up to r = 0.5a. However we expect that moving on 
the high density or low density regions (where the quartic polynomial becomes more steep) 
the quality of the fit will get worst. In particular the predicted negative values for an, 
in those regions, are completely unphysical and the fit should not be used to extrapolate 
beyond the range 0.28 < p < 0.65. 

IV. FROM WR TO NON ADDITIVE HARD SPHERES 

In order to see how the structure, and in particular the dcf's of the Widom-Rowlinson 
model change as one switches on the spheres diameters we have made two additional Monte 
Carlo simulations. In the first one we chose pi = p2 = O.Q5/Rf2 and -Rn = R22 = -R12/2. 
The resulting partial pair distribution functions and partial direct correlation functions are 
shown in Fig. ^ From a comparison with Fig. El we see how in this case the switching on of 
the like diameters causes both Ci2(r) for r < R12 and guir) for r > R12 to approach r = R12 
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with a slope close to zero. 

In the second simulation we chose pi = P2 = O.GS/i^fg -^u = O5 -^22 = Ri2- The 
resulting partial pair distribution functions and partial direct correlation functions are shown 
in Fig. 13 From a comparison with Fig. IHlwe see how in this case the switching on of the 
like diameters causes both (711(0) and cii(O) to increase, and cuir) to lose the nearly zero 
slope at r = 0. As in the previous case gi2{r) for r > R12 approaches r = R12 with a slope 
close to zero. The like 22 correlation functions for r > R12 vary over a range comparable to 
the one over which vary the like 11 correlation functions of the WR model. 

For both these cases there is no analytic solution of the PY approximation available and 
a better understanding of the behavior of the direct correlation functions may help in finding 
approximate expressions [l^ . 



V. CONCLUSIONS 



In this paper we have evaluated the direct correlation functions of a Widom-Rowlinson 
mixture at different densities through Monte Carlo simulation and we have studied the 
possibility of fitting the difference between MC data and the PY dcf 's. We found a very 
good parameterization of Ci2(r) for r < a [see equations (fT7j) and (fT8|l - |T9|) ] and a poorer one 
for cii(r) [see equations (|TT)|l and pn|) - P^ ]. The difficulty in this last case probably arises 
from the necessity of using three parameters [instead of just two needed for parameterizing 
Ci2(r)], although it cannot be completely excluded some effect of the decreasing precision of 
the simulation data near the origin. 

In the last part of the paper we have illustrated with additional Monte Carlo data the 
changes induced in the WR dcf 's by a finite size of the excluded volume of like correlations. 
These results are meant to provide a guide in the search of a manageable, simple analytical 
parameterization of the structure of mixtures of non additive hard spheres which is still not 
available although highly desirable. 
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LIST OF FIGURES 

Fig. ^ Top panel: partial direct correlation functions obtained from the Monte Carlo simula- 
tion (points) with the c"\2^\r) obtained from the PY approximation (line) at a density 
pa^ = 0.28748. Bottom panel: partial pair distribution functions obtained from the 
Monte Carlo simulation compared with the ones obtained from the PY approxima- 
tion at the same density. The open circles and the dashed line: the like correlation 
functions. Closed circles and the continuous line: the unlike correlation functions. 

Fig. 121 Same as in Fig.l at a density pa^ = 0.4. 

Fig. ini Same as in Fig.l at a density pa^ = 0.65. 

Fig. m We plot, for five different values of the density, the parameters ai2 (diagonal crosses) 
and bi2 (starred crosses) of the best Gaussian fit (fT7|) to A^2('") r < a, and fit 
them with parabolae (lines). The parameters at pa^ = 0.45 where not used for the 
parabolic fit and give an indication of the quality of the fit. 

Fig. El We plot, for five different values of the density, the parameters an, bu and dn (stars) of 
the best Gaussian fit |TH|) to A'j;]^(r), and draw the quartic polynomial (lines) through 
them. The parameters at pa^ = 0.45 where not used to determine the quartic poly- 
nomial and give an indication of the quality of the fit. 

Fig. ini Monte Carlo results at a density p = pi = P2 = 0.65/ for the partial direct 
correlation function (on top) and the partial pair distribution function (below) of a 
mixture of non additive hard spheres with Rn = R22 = -R12/2. The open circles 
denote the like correlation functions. The closed circles denote the unlike correlation 
functions. 

Fig. [71 Monte Carlo results at a density p = pi = P2 = 0.65/-Rf2 ^^e partial direct 
correlation function (on top) and the partial pair distribution function (below) of 
a mixture of non additive hard spheres with Ru = and R22 = Ri2- The open 
circles denote the like 11 correlation functions. The open triangles denote the like 22 
correlation functions. The closed circles denote the unlike correlation functions. 
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FIG. 1: R. Fantoni and G. Pastore 




FIG. 2: R. Fantoni and G. Pastore 
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FIG. 3: R. Fantoni and G. Pastore 
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FIG. 6: R. Fantoni and G. Pastore 




FIG. 7: R. Fantoni and G. Pastore 



